Translational Efficiency
###Comparison genotypes 14B and 4G, Contrast A
riboCountsA <- riboCounts %>% select(-contains("COL"))
rnaCountsA <- rnaCounts %>% select(-contains("COL"))
sampleAnnotationA <- sampleAnnotation %>% filter(!str_detect(group, 'COL'))
sampleAnnotation2A <- sampleAnnotation2 %>% filter(!str_detect(group, 'COL'))
# rna and ribo
sampleAnnotationA$SeqType = "RNA"
sampleAnnotation2A$SeqType = "Ribo"
combinedCountsA = cbind(riboCountsA, rnaCountsA)
sampleAnnotation3A = rbind(sampleAnnotationA, sampleAnnotation2A)
colnames(combinedCountsA) = rownames(sampleAnnotation3A)
# time + genotype + time:genotype + SeqType + SeqType:time + SeqType:genotype + SeqType:time:genotype
DESeqDataSet = DESeqDataSetFromMatrix(
countData = combinedCountsA,
colData = sampleAnnotation3A,
design = ~ time * genotype * SeqType
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet_both = DESeq(
DESeqDataSet,
parallel=FALSE,
test = "LRT",
reduced = ~ (time + genotype + SeqType)^2
)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results_both <- results(DESeqDataSet_both)
clean_DESeq_padj <- which(!is.na(DESeq_Results_both$padj))
sum(DESeq_Results_both[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 206
sum(DESeq_Results_both[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 20.6
# rna
DESeqDataSet = DESeqDataSetFromMatrix(
countData = rnaCountsA,
colData = sampleAnnotationA,
design = ~ time + genotype + time:genotype
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet = DESeq(
DESeqDataSet,
parallel=FALSE,
test = "LRT",
reduced = ~ time + genotype
)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results_RNA <- results(DESeqDataSet)
# ribo
DESeqDataSet = DESeqDataSetFromMatrix(
countData = riboCountsA,
colData = sampleAnnotation2A,
design = ~ time + genotype + time:genotype
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet = DESeq(
DESeqDataSet,
parallel=FALSE,
test = "LRT",
reduced = ~ time + genotype
)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results_ribo <- results(DESeqDataSet)
lgNorm = log2(counts(DESeqDataSet_both, normalized=TRUE) + 1)
Overall PCA Plot
pca = prcomp(t(lgNorm))
pcaData = data.frame(pca$x[ , 1:2])
pcaData$group = sampleAnnotation3A[rownames(pcaData), "group"]
pcaData$sample = rownames(pcaData)
gg = ggplot(pcaData, aes(x=PC1, y=PC2, color=group, label=sample))
gg = gg + geom_point(size=2.5, alpha=0.8)
gg = gg + scale_color_manual(values=groupColors)
gg = gg + theme(panel.background = element_rect(fill = 'aliceblue')) + ggtitle("Overall PCA for Translational Efficiency")
print(gg)

clean_DESeq_padj <- which(!is.na(DESeq_Results_both$padj))
clean_DESeq_Results <- DESeq_Results_both[clean_DESeq_padj, ]
significant_genes <- rownames(clean_DESeq_Results[which(clean_DESeq_Results[,"padj"] <= 0.1), ])
clean_DESeq_Results <- data.frame(clean_DESeq_Results[which(clean_DESeq_Results[,"padj"] <= 0.1), ])
clean_DESeq_Results$gene <- significant_genes
RNA_Result_groupings <- clean_DESeq_Results %>% inner_join(goAssociations2, by = "gene")
RNA_Result_groupings %>% group_by(gene_ontology_name) %>% summarize(count = n()) %>% arrange(desc(count)) %>%
filter(count >= 5) %>%
knitr::kable(caption = "Distribution of Significant Genes Groupings For TE")
Distribution of Significant Genes Groupings For TE
| response to light stimulus |
14 |
| translation |
14 |
| pyruvate metabolic process |
9 |
| response to cold |
6 |
| protein-chromophore linkage |
5 |
| translational elongation |
5 |
exclusive = rownames(DESeq_Results_both)[which(DESeq_Results_both$padj < 0.1 & DESeq_Results_ribo$padj < 0.1 & DESeq_Results_RNA$padj > 0.1)]
both = rownames(DESeq_Results_both)[which(DESeq_Results_both$padj < 0.1 & DESeq_Results_ribo$padj < 0.1 & DESeq_Results_RNA$padj < 0.1)]
intensified = both[which(DESeq_Results_both[both,2]*DESeq_Results_RNA[both,2] > 0)]
buffered = both[which(DESeq_Results_both[both,2]*DESeq_Results_RNA[both,2] < 0)]
goAssociations2 %>% filter(gene %in% exclusive) %>% group_by(gene_ontology_name) %>% summarise(count = n()) %>%
arrange(desc(count)) %>% knitr::kable(caption = "Gene Groupings for Exclusive TE Results")
Gene Groupings for Exclusive TE Results
| translation |
8 |
| response to light stimulus |
7 |
| NADH metabolic process |
2 |
| protein-chromophore linkage |
2 |
| proton motive force-driven ATP synthesis |
2 |
| response to cold |
2 |
| cold acclimation |
1 |
| glucose metabolic process |
1 |
| photosynthesis |
1 |
| pyruvate metabolic process |
1 |
| tetrapyrrole metabolic process |
1 |
goAssociations2 %>% filter(gene %in% both) %>% group_by(gene_ontology_name) %>% summarise(count = n()) %>% arrange(desc(count)) %>% knitr::kable(caption = "Gene Groupings for 'Both' TE Results")
Gene Groupings for ‘Both’ TE Results
| gluconeogenesis |
1 |
| pyruvate metabolic process |
1 |
| response to light stimulus |
1 |
goAssociations2 %>% filter(gene %in% intensified) %>% group_by(gene_ontology_name) %>% summarise(count = n()) %>% arrange(desc(count)) %>% knitr::kable(caption = "Gene Groupings for Intensified TE Results")
Gene Groupings for Intensified TE Results
| gluconeogenesis |
1 |
| pyruvate metabolic process |
1 |
goAssociations2 %>% filter(gene %in% buffered) %>% group_by(gene_ontology_name) %>% summarise(count = n()) %>% arrange(desc(count)) %>% knitr::kable(caption = "Gene Groupings for Buffered TE Results")
Gene Groupings for Buffered TE Results
| response to light stimulus |
1 |
goAssociations2 %>% filter(gene %in% exclusive) %>%
knitr::kable(caption = "Gene Groupings for Exclusive TE Results")
goAssociations2 %>% filter(gene %in% both) %>%
knitr::kable(caption = "Gene Groupings for 'Both' TE Results")
Gene Groupings for ‘Both’ TE Results
| GO:0006090 |
pyruvate metabolic process |
AT3G12780 |
| GO:0006094 |
gluconeogenesis |
AT3G12780 |
| GO:0009416 |
response to light stimulus |
AT1G51200 |
goAssociations2 %>% filter(gene %in% intensified) %>%
knitr::kable(caption = "Gene Groupings for Intensified TE Results")
Gene Groupings for Intensified TE Results
| GO:0006090 |
pyruvate metabolic process |
AT3G12780 |
| GO:0006094 |
gluconeogenesis |
AT3G12780 |
goAssociations2 %>% filter(gene %in% buffered) %>%
knitr::kable(caption = "Gene Groupings for Buffered TE Results")
Gene Groupings for Buffered TE Results
| GO:0009416 |
response to light stimulus |
AT1G51200 |
for (id in both){
ymax=max(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
ymin=min(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
plot(c(0,1), c(0,DESeq_Results_ribo[id,2]), type="l",col="gray", ylim=c(ymin,ymax),
ylab="Log2 fold change",xlab="",xaxt="n")
lines(c(0,1), c(0,DESeq_Results_RNA[id,2]),type="l",col="blue")
lines(c(0,1), c(0,DESeq_Results_both[id,2]), type="l",col="red")
legend("bottomleft",c("RNA","Ribo","TE"),fill=c("blue","gray","red"), cex=1, border = NA, bty="n")
axis(1,at=c(0,1),labels=c(1,2),las=1)
title(id)
}

Intensified and buffered: Genes regulated both by transcriptional and
translational regulation (significant ΔRNA, ΔRPFs, and ΔTE) include
intensified and buffered genes. These genes are both DTGs and DTEGs.
All lines going in the same direction –> change in translational
efficiency is counteracting the change in RNA
for (id in exclusive){
ymax=max(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
ymin=min(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
plot(c(0,1), c(0,DESeq_Results_ribo[id,2]), type="l",col="gray", ylim=c(ymin,ymax),
ylab="Log2 fold change",xlab="",xaxt="n")
lines(c(0,1), c(0,DESeq_Results_RNA[id,2]),type="l",col="blue")
lines(c(0,1), c(0,DESeq_Results_both[id,2]), type="l",col="red")
legend("bottomleft",c("RNA","Ribo","TE"),fill=c("blue","gray","red"), cex=1, border = NA, bty="n")
axis(1,at=c(0,1),labels=c(1,2),las=1)
title(id)
}



























exclusive focuses on findings that are translationally different
only.
for (id in intensified){
ymax=max(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
ymin=min(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
plot(c(0,1), c(0,DESeq_Results_ribo[id,2]), type="l",col="gray", ylim=c(ymin,ymax),
ylab="Log2 fold change",xlab="",xaxt="n")
lines(c(0,1), c(0,DESeq_Results_RNA[id,2]),type="l",col="blue")
lines(c(0,1), c(0,DESeq_Results_both[id,2]), type="l",col="red")
legend("bottomleft",c("RNA","Ribo","TE"),fill=c("blue","gray","red"), cex=1, border = NA, bty="n")
axis(1,at=c(0,1),labels=c(1,2),las=1)
title(id)
}

for (id in buffered){
ymax=max(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
ymin=min(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
plot(c(0,1), c(0,DESeq_Results_ribo[id,2]), type="l",col="gray", ylim=c(ymin,ymax),
ylab="Log2 fold change",xlab="",xaxt="n")
lines(c(0,1), c(0,DESeq_Results_RNA[id,2]),type="l",col="blue")
lines(c(0,1), c(0,DESeq_Results_both[id,2]), type="l",col="red")
legend("bottomleft",c("RNA","Ribo","TE"),fill=c("blue","gray","red"), cex=1, border = NA, bty="n")
axis(1,at=c(0,1),labels=c(1,2),las=1)
title(id)
}

###Comparison genotypes Col and 14B, Contrast B
riboCountsB <- riboCounts %>% select(-contains("4G"))
rnaCountsB <- rnaCounts %>% select(-contains("4G"))
sampleAnnotationB <- sampleAnnotation %>% filter(!str_detect(group, '4G'))
sampleAnnotation2B <- sampleAnnotation2 %>% filter(!str_detect(group, '4G'))
# rna and ribo
sampleAnnotationB$SeqType = "RNA"
sampleAnnotation2B$SeqType = "Ribo"
combinedCountsB = cbind(riboCountsB, rnaCountsB)
sampleAnnotation3B = rbind(sampleAnnotationB, sampleAnnotation2B)
colnames(combinedCountsB) = rownames(sampleAnnotation3B)
# time + genotype + time:genotype + SeqType + SeqType:time + SeqType:genotype + SeqType:time:genotype
DESeqDataSet = DESeqDataSetFromMatrix(
countData = combinedCountsB,
colData = sampleAnnotation3B,
design = ~ time * genotype * SeqType
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet_both = DESeq(
DESeqDataSet,
parallel=FALSE,
test = "LRT",
reduced = ~ (time + genotype + SeqType)^2
)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results_both <- results(DESeqDataSet_both)
clean_DESeq_padj <- which(!is.na(DESeq_Results_both$padj))
sum(DESeq_Results_both[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 145
sum(DESeq_Results_both[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 14.5
# rna
DESeqDataSet = DESeqDataSetFromMatrix(
countData = rnaCountsB,
colData = sampleAnnotationB,
design = ~ time + genotype + time:genotype
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet = DESeq(
DESeqDataSet,
parallel=FALSE,
test = "LRT",
reduced = ~ time + genotype
)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results_RNA <- results(DESeqDataSet)
# ribo
DESeqDataSet = DESeqDataSetFromMatrix(
countData = riboCountsB,
colData = sampleAnnotation2B,
design = ~ time + genotype + time:genotype
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet = DESeq(
DESeqDataSet,
parallel=FALSE,
test = "LRT",
reduced = ~ time + genotype
)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results_ribo <- results(DESeqDataSet)
lgNorm = log2(counts(DESeqDataSet_both, normalized=TRUE) + 1)
Overall PCA Plot
pca = prcomp(t(lgNorm))
pcaData = data.frame(pca$x[ , 1:2])
pcaData$group = sampleAnnotation3B[rownames(pcaData), "group"]
pcaData$sample = rownames(pcaData)
gg = ggplot(pcaData, aes(x=PC1, y=PC2, color=group, label=sample))
gg = gg + geom_point(size=2.5, alpha=0.8)
gg = gg + scale_color_manual(values=groupColors)
gg = gg + theme(panel.background = element_rect(fill = 'aliceblue')) + ggtitle("Overall PCA for Translational Efficiency")
print(gg)

clean_DESeq_padj <- which(!is.na(DESeq_Results_both$padj))
clean_DESeq_Results <- DESeq_Results_both[clean_DESeq_padj, ]
significant_genes <- rownames(clean_DESeq_Results[which(clean_DESeq_Results[,"padj"] <= 0.1), ])
clean_DESeq_Results <- data.frame(clean_DESeq_Results[which(clean_DESeq_Results[,"padj"] <= 0.1), ])
clean_DESeq_Results$gene <- significant_genes
RNA_Result_groupings <- clean_DESeq_Results %>% inner_join(goAssociations2, by = "gene")
RNA_Result_groupings %>% group_by(gene_ontology_name) %>% summarize(count = n()) %>% arrange(desc(count)) %>%
filter(count >= 5) %>%
knitr::kable(caption = "Distribution of Significant Genes Groupings For TE")
Distribution of Significant Genes Groupings For TE
| pyruvate metabolic process |
9 |
| translation |
9 |
| photosynthesis |
8 |
| response to light stimulus |
8 |
| protein-chromophore linkage |
6 |
| response to cold |
6 |
exclusive = rownames(DESeq_Results_both)[which(DESeq_Results_both$padj < 0.1 & DESeq_Results_ribo$padj < 0.1 & DESeq_Results_RNA$padj > 0.1)]
both = rownames(DESeq_Results_both)[which(DESeq_Results_both$padj < 0.1 & DESeq_Results_ribo$padj < 0.1 & DESeq_Results_RNA$padj < 0.1)]
intensified = both[which(DESeq_Results_both[both,2]*DESeq_Results_RNA[both,2] > 0)]
buffered = both[which(DESeq_Results_both[both,2]*DESeq_Results_RNA[both,2] < 0)]
goAssociations2 %>% filter(gene %in% exclusive) %>% group_by(gene_ontology_name) %>% summarise(count = n()) %>%
arrange(desc(count)) %>% knitr::kable(caption = "Gene Groupings for Exclusive TE Results")
Gene Groupings for Exclusive TE Results
| photosynthesis |
7 |
| protein-chromophore linkage |
6 |
| response to light stimulus |
5 |
| translation |
4 |
| response to cold |
3 |
| photosynthetic electron transport chain |
2 |
| proton motive force-driven ATP synthesis |
2 |
| transmembrane transport |
2 |
| cold acclimation |
1 |
| glucose metabolic process |
1 |
| negative regulation of catalytic activity |
1 |
| protein ubiquitination |
1 |
| pyruvate metabolic process |
1 |
| response to disaccharide |
1 |
| response to heat |
1 |
| response to ozone |
1 |
| tetrapyrrole metabolic process |
1 |
| thylakoid membrane organization |
1 |
| valine biosynthetic process |
1 |
goAssociations2 %>% filter(gene %in% both) %>% group_by(gene_ontology_name) %>% summarise(count = n()) %>% arrange(desc(count)) %>% knitr::kable(caption = "Gene Groupings for 'Both' TE Results")
Gene Groupings for ‘Both’ TE Results
| translation |
3 |
| translational elongation |
2 |
| gluconeogenesis |
1 |
| pyruvate metabolic process |
1 |
goAssociations2 %>% filter(gene %in% intensified) %>% group_by(gene_ontology_name) %>% summarise(count = n()) %>% arrange(desc(count)) %>% knitr::kable(caption = "Gene Groupings for Intensified TE Results")
Gene Groupings for Intensified TE Results
| gluconeogenesis |
1 |
| pyruvate metabolic process |
1 |
| translation |
1 |
goAssociations2 %>% filter(gene %in% buffered) %>% group_by(gene_ontology_name) %>% summarise(count = n()) %>% arrange(desc(count)) %>% knitr::kable(caption = "Gene Groupings for Buffered TE Results")
Gene Groupings for Buffered TE Results
| translation |
2 |
| translational elongation |
2 |
goAssociations2 %>% filter(gene %in% exclusive) %>%
knitr::kable(caption = "Gene Groupings for Exclusive TE Results")
goAssociations2 %>% filter(gene %in% both) %>%
knitr::kable(caption = "Gene Groupings for 'Both' TE Results")
goAssociations2 %>% filter(gene %in% intensified) %>%
knitr::kable(caption = "Gene Groupings for Intensified TE Results")
Gene Groupings for Intensified TE Results
| GO:0006090 |
pyruvate metabolic process |
AT3G12780 |
| GO:0006094 |
gluconeogenesis |
AT3G12780 |
| GO:0006412 |
translation |
AT3G27160 |
goAssociations2 %>% filter(gene %in% buffered) %>%
knitr::kable(caption = "Gene Groupings for Buffered TE Results")
for (id in both){
ymax=max(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
ymin=min(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
plot(c(0,1), c(0,DESeq_Results_ribo[id,2]), type="l",col="gray", ylim=c(ymin,ymax),
ylab="Log2 fold change",xlab="",xaxt="n")
lines(c(0,1), c(0,DESeq_Results_RNA[id,2]),type="l",col="blue")
lines(c(0,1), c(0,DESeq_Results_both[id,2]), type="l",col="red")
legend("bottomleft",c("RNA","Ribo","TE"),fill=c("blue","gray","red"), cex=1, border = NA, bty="n")
axis(1,at=c(0,1),labels=c(1,2),las=1)
title(id)
}




Intensified and buffered: Genes regulated both by transcriptional and
translational regulation (significant ΔRNA, ΔRPFs, and ΔTE) include
intensified and buffered genes. These genes are both DTGs and DTEGs.
All lines going in the same direction –> change in translational
efficiency is counteracting the change in RNA
for (id in exclusive){
ymax=max(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
ymin=min(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
plot(c(0,1), c(0,DESeq_Results_ribo[id,2]), type="l",col="gray", ylim=c(ymin,ymax),
ylab="Log2 fold change",xlab="",xaxt="n")
lines(c(0,1), c(0,DESeq_Results_RNA[id,2]),type="l",col="blue")
lines(c(0,1), c(0,DESeq_Results_both[id,2]), type="l",col="red")
legend("bottomleft",c("RNA","Ribo","TE"),fill=c("blue","gray","red"), cex=1, border = NA, bty="n")
axis(1,at=c(0,1),labels=c(1,2),las=1)
title(id)
}





































exclusive focuses on findings that are translationally different
only.
for (id in intensified){
ymax=max(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
ymin=min(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
plot(c(0,1), c(0,DESeq_Results_ribo[id,2]), type="l",col="gray", ylim=c(ymin,ymax),
ylab="Log2 fold change",xlab="",xaxt="n")
lines(c(0,1), c(0,DESeq_Results_RNA[id,2]),type="l",col="blue")
lines(c(0,1), c(0,DESeq_Results_both[id,2]), type="l",col="red")
legend("bottomleft",c("RNA","Ribo","TE"),fill=c("blue","gray","red"), cex=1, border = NA, bty="n")
axis(1,at=c(0,1),labels=c(1,2),las=1)
title(id)
}


for (id in buffered){
ymax=max(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
ymin=min(DESeq_Results_ribo[id,2],DESeq_Results_RNA[id,2],DESeq_Results_both[id,2],0)
plot(c(0,1), c(0,DESeq_Results_ribo[id,2]), type="l",col="gray", ylim=c(ymin,ymax),
ylab="Log2 fold change",xlab="",xaxt="n")
lines(c(0,1), c(0,DESeq_Results_RNA[id,2]),type="l",col="blue")
lines(c(0,1), c(0,DESeq_Results_both[id,2]), type="l",col="red")
legend("bottomleft",c("RNA","Ribo","TE"),fill=c("blue","gray","red"), cex=1, border = NA, bty="n")
axis(1,at=c(0,1),labels=c(1,2),las=1)
title(id)
}



###Comparison genotypes Col and 4G, Contrast C
riboCountsC <- riboCounts %>% select(-contains("14B"))
rnaCountsC <- rnaCounts %>% select(-contains("14B"))
sampleAnnotationC <- sampleAnnotation %>% filter(!str_detect(group, '14B'))
sampleAnnotation2C <- sampleAnnotation2 %>% filter(!str_detect(group, '14B'))
# rna and ribo
sampleAnnotationC$SeqType = "RNA"
sampleAnnotation2C$SeqType = "Ribo"
combinedCountsC = cbind(riboCountsC, rnaCountsC)
sampleAnnotation3C = rbind(sampleAnnotationC, sampleAnnotation2C)
colnames(combinedCountsC) = rownames(sampleAnnotation3C)
# time + genotype + time:genotype + SeqType + SeqType:time + SeqType:genotype + SeqType:time:genotype
DESeqDataSet = DESeqDataSetFromMatrix(
countData = combinedCountsC,
colData = sampleAnnotation3C,
design = ~ time * genotype * SeqType
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet_both = DESeq(
DESeqDataSet,
parallel=FALSE,
test = "LRT",
reduced = ~ (time + genotype + SeqType)^2
)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results_both <- results(DESeqDataSet_both)
clean_DESeq_padj <- which(!is.na(DESeq_Results_both$padj))
sum(DESeq_Results_both[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 0
sum(DESeq_Results_both[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 0
# rna
DESeqDataSet = DESeqDataSetFromMatrix(
countData = rnaCountsC,
colData = sampleAnnotationC,
design = ~ time + genotype + time:genotype
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet = DESeq(
DESeqDataSet,
parallel=FALSE,
test = "LRT",
reduced = ~ time + genotype
)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results_RNA <- results(DESeqDataSet)
# ribo
DESeqDataSet = DESeqDataSetFromMatrix(
countData = riboCountsC,
colData = sampleAnnotation2C,
design = ~ time + genotype + time:genotype
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet = DESeq(
DESeqDataSet,
parallel=FALSE,
test = "LRT",
reduced = ~ time + genotype
)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results_ribo <- results(DESeqDataSet)
lgNorm = log2(counts(DESeqDataSet_both, normalized=TRUE) + 1)
Overall PCA Plot
pca = prcomp(t(lgNorm))
pcaData = data.frame(pca$x[ , 1:2])
pcaData$group = sampleAnnotation3C[rownames(pcaData), "group"]
pcaData$sample = rownames(pcaData)
gg = ggplot(pcaData, aes(x=PC1, y=PC2, color=group, label=sample))
gg = gg + geom_point(size=2.5, alpha=0.8)
gg = gg + scale_color_manual(values=groupColors)
gg = gg + theme(panel.background = element_rect(fill = 'aliceblue')) + ggtitle("Overall PCA for Translational Efficiency")
print(gg)
